%% Calculate equilibrium

% function [P_hat_mnj,w_hat,X_hat_mnj,pi_mnjf1,pi_mnj1,...
%     pi_c_mnj1,pi_c_nj1,square_brackets,P_hat_n,...
%     pseudo_final_goods,pseudo_intermediates,pseudo_intermediates_mnjj,...
%     pi_M1,pi_l1,pi_M_f1,pi_l_f1]= ...
%     EquilibriumCES_fun(P_hat_mnj,w_hat,pi_mnj0,pi_M0,pi_M_f0,...
%     pi_l0,pi_l_f0,pi_mnjf0,pi_c_mnj0,X_hat_mnj,pi_c_nj0,PC_n,zeta_mnj)

function [P_hat_n,P_hat_nj,P_hat_mnj,w_hat,X_hat_mnj,pi_mnjf1,pi_mnj1,...
    pi_c_mnj1,pi_c_nj1, pi_M1,pi_l1,pi_M_f1,pi_l_f1,PC_hat_n1,...
    pseudo_final_goods1,pseudo_intermediates1,mu_hat_mnjf]= ...
    EquilibriumCES_fun_approx_Olig(P_hat_mnj0,w_hat,pi_mnj0,pi_M0,pi_M_f0,...
    pi_l0,pi_l_f0,pi_mnjf0,pi_c_mnj0,X_hat_mnj0,pi_c_nj0,PC_n,zeta_mnj,X_mnj,mu_hat_mnjf0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF  
%
% Program: Guess prices and wages, then iterate over prices, compare this 
% with the market clearing condition for the new guess of wages, then 
% iterate until convergence.
%
% Calls up the following functions:
%  1) PriceLoopCES_fun.m
% 2) TradeShareCES_fun.m
% 	3) ExpenditureCES_fun.m
%	4) LabourMarketConditionCES_fun.m
% 
% Output: Converged expenditure hat, wage hat and prices 
% 
% NEW: For all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}. Global variables 
% changed and slight modification to file structure for calling up data.
% We also get rid of gamma_type variable. What previously _2.csv. files
% slightly as well as adjusting in our new directory.
%
% We also get rid of using the .(ISO{}) cell call up for the gammas and
% instead use a large matrix saved in Step2. This will help save on loops
% and speed up the proceedure.
%
% Last Updated: 5/02/2019 by IM & JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma lambda eta ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej ...
    sector_algorithm 

w_max = 1; e = 1;

W = [w_hat'];
D = [];

fprintf('WAGE ITERATION=%d \n',e);

while (e <= maxit) && (w_max > tol)
   % Price Loop: Takes wage and price guess and loops until price
   % converges. Output is is the new P hat (change in prices) and wage hat
   % as well as expenditure hat (change in inputs/exports).
      
   [P_hat_mnj,b_hat,b_hat_f,P_M_hat,P_M_f_hat,mu_hat_mnjf] = PriceLoopCES_fun_approx_Olig(P_hat_mnj0,w_hat,pi_mnj0,...
          pi_M0,pi_l0,pi_mnjf0,pi_M_f0,pi_l_f0,mu_hat_mnjf0);

   % Trade and labor Share: Creates the new trade (export and intermediate
   % input) and labor shares for the next period (t+1). Note that 
   % intermediate input and labor shares are not fed into P_hat_mnj 
   % function abovve, which is still a function of initial values only

 [pi_mnj1,pi_mnjf1,pi_c_mnj1,pi_M1,pi_l1,pi_M_f1,pi_l_f1] = ...
   TradeShareCESOlig_fun(P_hat_mnj,b_hat,b_hat_f,P_M_f_hat,pi_M0,pi_l0,w_hat,pi_mnj0,...
   pi_M_f0,pi_l_f0,pi_mnjf0,pi_c_mnj0,mu_hat_mnjf);

   % Expenditure: Creates a guess for the change in expenditure from t to 
   % t+1 and then loops over this until it converges. Follows the 
   % expenditure hat equation 23 from the master notes
   
   [P_hat_n,P_hat_nj,X_hat_mnj,pi_c_nj1,PC_hat_n1,...
       pseudo_intermediates1,pseudo_final_goods1] = ...
      ExpenditureCESOlig_fun(P_hat_mnj,pi_M1,pi_l1,w_hat,...
     pi_M_f1,pi_l_f1,pi_c_mnj0,PC_n,pi_c_nj0,...
   X_mnj,pi_c_mnj1,zeta_mnj,pi_mnjf1,pi_mnj1,X_hat_mnj0);
  
   % Labour Market Condition: Uses the labour market condition to update
   %                          the wage guess
      
   [w_hat_new] = LabourMarketConditionCESOlig_fun(pi_l_f0, ...
            pi_l0,pi_l_f1,pi_l1,X_mnj,X_hat_mnj,...
            P_hat_n,pi_mnjf0,pi_mnjf1,pi_mnj0,pi_mnj1);  
%  
%         
    P_hat_mnj0 = P_hat_mnj;
    X_hat_mnj0 = X_hat_mnj;
%    
%    %w_max=sum(abs(w_hat_new-w_hat));
%    
    w_diff=abs(w_hat_new-w_hat)./w_hat_new;
    w_max=norm(w_diff);
% 
%    
    if alphaT == 2
%        
%      % Convex iteration guess updater should be more stable and more likely
%      % to converge as smaller jumps in the guess, but is slower to converge.
%    
      w_hat=0.5*w_hat_new+0.5*w_hat;
%      
    elseif alphaT == 4
%        
%      % Convex iteration guess updater should be more stable and more likely
%      % to converge as smaller jumps in the guess, but is slower to converge.
%    
%     % w_hat=0.5*w_hat_new+0.5*w_hat;
      w_hat=w_hat_new; 
%      % Updating the wage guess, this is the fastest but a convex iteration
%      % updater may be more stable (smaller jumps of the guess)
    else
%        
     w_hat=w_hat_new;
%    
%     
    end
   
%    % Update shares
%    pi_mnj0 = pi_mnj1;
%    pi_mnjf0 = pi_mnjf1;
%    pi_M0 = pi_M1;
%    pi_M_f0 = pi_M_f1;
%    pi_l0 = pi_l1;
%    pi_l_f0 = pi_l_f1;
%    pi_c_mnj0 = pi_c_mnj1;
    
   
%    % Small interation wage relative to previous step given a nu_W, where
%    % nu_W is the adjustment we might make
%    
%    nu_W = .5;
%    
%    w_diff = (w_hat_new-w_hat)./w_hat;
%    
%    w_hat_new = w_hat + nu_W*w_hat.*w_diff; % Iteration of new wage
%    
%    w_diff1 = (w_hat_new-w_hat)./w_hat;
%    
%    % Calculate distance as max between non-adjusted and adjusted new w_hat 
%    % and old w_hawt
% 
%    w_max = max(max(abs(w_diff)),max(abs(w_diff1)));
   
   % Set wage to new value to feed into next loop iteration
   
     W = [W w_hat_new']
     D = [D w_max]
%    
% %   w_hat = w_hat_new;
%    
    e=e+1;
% 
    fprintf('WAGE ITERATION=%d        Max wage distance=%d  \n',e-1,w_max);
end

end